### Analysis Code for "Value Similarity and Norm Change: Null Effects and Backlash to Messaging on Same-Sex Rights in Uganda"
### Purpose: to produce analysis for main text results (except fig. 1)
### Author: Nicholas Lyon

# Preliminaries -----------------------------------------------------------

# Purpose of code section: load required libraries and cleaned data

library('rio')
library('tidyverse')
library('scales')
library('ggplot2')
library('ggthemes')
library('estimatr')
library('texreg')
library('list')

root <- "add your path here/" # add path to where "uganda_samesex_replication" is stored
setwd(paste0(root, "uganda_samesex_replication/Output"))

survey <- import(paste0(root, 'uganda_samesex_replication/Data/clean_data.csv')) # bring in cleaned survey data
list <- import(paste0(root, 'uganda_samesex_replication/Data/list_data.csv')) # bring in cleaned list analysis data


# Figure 2 ----------------------------------------------------------------

# Panel (a): value similarity with other African countries
ggplot(survey, aes(x = factor(values_africa))) + 
  geom_bar(aes(y = (..count..)/sum(..count..)), fill = 'black') +
  geom_text(aes(label = paste0(100*round((..count..)/sum(..count..), 
                                         2), "%"), 
                y= (..count..)/sum(..count..)), 
            stat= "count", vjust = -.5) +
  scale_y_continuous(labels = percent) + 
  scale_x_discrete(labels=c("0" = "Not at all",
                            "1" = "Not very",
                            "2" = "Somewhat",
                            "3" = "Very",
                            "4" = "Same",
                            "99" = "Not \nsure")) + 
  coord_cartesian(ylim = c(0, .7)) +
  labs(x = "\nHow Similar?",
       y = "Percent \n") + 
  theme_clean() + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1),
        axis.text=element_text(size=12),
        axis.title=element_text(size=14, face="bold"),
        legend.title = element_blank()) +
  ggsave(filename = "african_plot.pdf", 
         device = "pdf", width = 4, height = 3, units = "in")

# Panel (b): value similarity with Western countries
ggplot(survey, aes(x = factor(values_west))) + 
  geom_bar(aes(y = (..count..)/sum(..count..)), fill = 'black') +
  geom_text(aes(label = paste0(100*round((..count..)/sum(..count..), 
                                         2), "%"), 
                y= (..count..)/sum(..count..)), 
            stat= "count", vjust = -.5) +
  scale_y_continuous(labels = percent) + 
  scale_x_discrete(labels=c("0" = "Not at all",
                            "1" = "Not very",
                            "2" = "Somewhat",
                            "3" = "Very",
                            "4" = "Same",
                            "99" = "Not \nsure")) + 
  coord_cartesian(ylim = c(0, .7)) +
  labs(x = "\nHow Similar?",
       y = "Percent \n") + 
  theme_clean() + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1),
        axis.text=element_text(size=12),
        axis.title=element_text(size=14, face="bold"),
        legend.title = element_blank()) +
  ggsave(filename = "western_plot.pdf", 
         device = "pdf", width = 4, height = 3, units = "in")

# Panel (c): value similarity comparison
ggplot(subset(survey, !is.na(survey$values_compare)), aes(x = factor(values_compare))) + 
  # above line filters off 6 missing values from refused to answer responses
  geom_bar(aes(y = (..count..)/sum(..count..)), fill = 'black') +
  geom_text(aes(label = paste0(100*round((..count..)/sum(..count..), 
                                         2), "%"), 
                y= (..count..)/sum(..count..)), 
            stat= "count", vjust = -.5) +
  scale_y_continuous(labels = percent) + 
  scale_x_discrete(labels=c("1" = "Other \nAfrican \ncountries",
                            "0" = "About \nthe \nsame",
                            "-1" = "Non-African \ncountries",
                            "99" = "Not \nsure")) + 
  labs(x = "\nWho Has More Similar Values?",
       y = "Percent \n") + 
  coord_cartesian(ylim = c(0, .85)) +
  theme_clean() + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1),
        axis.text=element_text(size=12),
        axis.title=element_text(size=14, face="bold"),
        legend.title = element_blank()) +
  ggsave(filename = "compare_plot.pdf", 
         device = "pdf", width = 4, height = 3, units = "in")

# Table 2 -----------------------------------------------------------------

# relevel treatment such that placebo is reference category for regressions
survey$treatment <- relevel(factor(survey$treatment), ref = 2)

# african: legal
ap_af_legal <- lm_robust(legal_africa ~ treatment, data = survey, se_type = "HC2")
# non-african:legal
ap_non_legal <- lm_robust(legal_nonafrican ~ treatment, data = survey, se_type = "HC2")
# african: welcoming
ap_af_attitude <- lm_robust(african_welcome_homosexuals ~ treatment, data = survey, se_type = "HC2")
# non-african: welcoming
ap_non_attitude <- lm_robust(nonafrican_welcome_homosexuals ~ treatment, data = survey, se_type = "HC2")

# combine into table
texreg(list(ap_af_legal, ap_non_legal, ap_af_attitude, ap_non_attitude), include.ci = FALSE,
       custom.coef.map = list("(Intercept)" = "Intercept",
                              "treatmentafrican" = "African",
                              "treatmentwestern" = "Western"),
       custom.model.names = c("African: Legal", "Non-African: Legal", 
                              "African: Welcoming", "Non-African: Welcoming"),
       digits = 3,
       include.rmse = FALSE,
       caption = "Manipulation Checks",
       stars = c(0.001, 0.01, 0.05))

# Table 3 -----------------------------------------------------------------

# social distance
reg_sd <- lm_robust(social_distance ~ treatment, data = survey, se_type = "HC2")
# rights index
reg_ss_rights_index <- lm_robust(ss_rights_index ~ treatment, data = survey, se_type = "HC2")
# morality index
reg_ss_morality_index <- lm_robust(ss_morality_index ~ treatment, data = survey, se_type = "HC2")
# mp message
reg_message <- lm_robust(homosexual_message ~ treatment, data = survey, se_type = "HC2")
# contribution game
reg_destruction_game <- lm_robust(destruction_game ~ treatment, data = survey, se_type = "HC2")


# combine into table
texreg(list(reg_sd, reg_ss_rights_index, reg_ss_morality_index, reg_message, reg_destruction_game), 
       include.ci = FALSE,
       custom.coef.map = list("(Intercept)" = "Intercept",
                              "treatmentafrican" = "African",
                              "treatmentwestern" = "Western"),
       custom.model.names = c("Social Distance", "Rights Index", 
                              "Morality Index", "MP Message", "Contribution Game"),
       digits = 3,
       include.rmse = FALSE,
       caption = "Treatment Effects of Messaging",
       stars = c(0.001, 0.01, 0.05))

# Table 4 -----------------------------------------------------------------

# using ictreg from "list" package (Blair and Imai) for each individual list estimate
# the estimates from each list are combined using simple average of estimates from both lists
# and standard errors on the double list estimate are produced using the method in Droitcour et al., defined below

droitcour_se <- function(list1, list1_t, list2, list2_t){
  x4a <- list1[list1_t == 1]
  y3a <- list2[list2_t == 0]
  y4b <- list2[list2_t == 1]
  x3b <- list1[list1_t == 0]
  
  n <- length(x4a) + length(x3b)
  
  est_var <- 1/4*(var(x4a - y3a) + var(y4b - x3b))
  est_se <- sqrt(est_var) / sqrt(n)
  return(est_se)
}

# double list estimate and SE for placebo control
placebo <- subset(list, treatment == "placebo") # subset to placebo arm

placebo_list1 <- ictreg(list1 ~ 1, data = placebo, 
                        treat = "list1_t", J = 4, method = "lm") # estimate from list 1
placebo_list2 <- ictreg(list2 ~ 1, data = placebo, 
                        treat = "list2_t", J = 4, method = "lm") # estimate from list 2

placebo_estimate <- (placebo_list1$par.treat + placebo_list2$par.treat) / 2 # double list estimate
placebo_se <- droitcour_se(list1 = placebo$list1, list1_t = placebo$list1_t, 
                           list2 = placebo$list2, list2_t = placebo$list2_t) # double list se
placebo_p <- 2*pnorm(placebo_estimate / placebo_se, lower.tail = F) # double list p-value 

# double list estimate and SE for african message 
african <- subset(list, treatment == "african") # subset to african message arm

african_list1 <- ictreg(list1 ~ 1, data = african, 
                        treat = "list1_t", J = 4, method = "lm") # estimate from list 1
african_list2 <- ictreg(list2 ~ 1, data = african, 
                        treat = "list2_t", J = 4, method = "lm") # estimate from list 2

african_estimate <- (african_list1$par.treat + african_list2$par.treat) / 2 # double list estimate
african_se <- droitcour_se(list1 = african$list1, list1_t = african$list1_t, 
                           list2 = african$list2, list2_t = african$list2_t) # double list se
african_p <- 2*pnorm(african_estimate / african_se, lower.tail = F) # double list p-value 

# double list estimate and SE for western message 
western <- subset(list, treatment == "western") # subset to western message arm

western_list1 <- ictreg(list1 ~ 1, data = western, 
                        treat = "list1_t", J = 4, method = "lm") # estimate from list 1
western_list2 <- ictreg(list2 ~ 1, data = western, 
                        treat = "list2_t", J = 4, method = "lm") # estimate from list 2

western_estimate <- (western_list1$par.treat + western_list2$par.treat) / 2 # double list estimate
western_se <- droitcour_se(list1 = western$list1, list1_t = western$list1_t, 
                           list2 = western$list2, list2_t = western$list2_t) # double list se
western_p <- 2*pnorm(western_estimate / western_se, lower.tail = F) # double list p-value 

# difference in means: african message and control
diff_african_placebo_estimate <- african_estimate - placebo_estimate # diff-in-means estimate
diff_african_placebo_se <- sqrt(african_se^2 + placebo_se^2) # diff-in-means se
diff_african_placebo_p <- 2*pnorm(diff_african_placebo_estimate / diff_african_placebo_se, lower.tail = F) # diff-in-means p-value

# difference in means: western message and control
diff_western_placebo_estimate <- western_estimate - placebo_estimate # diff-in-means estimate
diff_western_placebo_se <- sqrt(western_se^2 + placebo_se^2) # diff-in-means se
diff_western_placebo_p <- 2*pnorm(diff_western_placebo_estimate / diff_western_placebo_se, lower.tail = T) # diff-in-means p-value

# difference in means: african message and western message
diff_african_western_estimate <- african_estimate - western_estimate # diff-in-means estimate
diff_african_western_se <- sqrt(african_se^2 + western_se^2) # diff-in-means se
diff_african_western_p <- 2*pnorm(diff_african_western_estimate / diff_african_western_se, lower.tail = F) # diff-in-means p-value


# Table 5 -----------------------------------------------------------------

# african country perceptions
african_countries <- lm_robust(african_perceptions ~ treatment, 
                              data = survey, se_type = "HC2")
# western country perceptions
western_countries <- lm_robust(nonafrican_perceptions ~ treatment, 
                             data = survey, se_type = "HC2")

# combine into table
texreg(list(african_countries, western_countries), 
       include.ci = FALSE,
       custom.coef.map = list("(Intercept)" = "Intercept",
                              "treatmentafrican" = "African",
                              "treatmentwestern" = "Western"),
       custom.model.names = c("African Country Perceptions", 
                              "Non-African Country Perceptions"),
       digits = 3,
       include.rmse = FALSE,
       caption = "Treatment Effects on Country Perceptions",
       stars = c(0.001, 0.01, 0.05))


# Figure 3 ----------------------------------------------------------------

# combine messaging effects into one visualization - join results into data frame
main_viz <- tibble(
  outcome = c("Social Distance", "Social Distance",
              "Rights Index", "Rights Index",
              "Morality Index", "Morality Index",
              "MP Message", "MP Message",
              "Contribution Game", "Contribution Game",
              "Double List", "Double List",
              "African Country Perceptions", "African Country Perceptions",
              "Western Country Perceptions", "Western Country Perceptions"),
  type = c("African Message", "Western Message",
           "African Message", "Western Message",
           "African Message", "Western Message",
           "African Message", "Western Message",
           "African Message", "Western Message",
           "African Message", "Western Message",
           "African Message", "Western Message",
           "African Message", "Western Message"),
  estimate = c(reg_sd$coefficients[2], reg_sd$coefficients[3],
               reg_ss_rights_index$coefficients[2], reg_ss_rights_index$coefficients[3],
               reg_ss_morality_index$coefficients[2], reg_ss_morality_index$coefficients[3],
               reg_message$coefficients[2], reg_message$coefficients[3],
               reg_destruction_game$coefficients[2], reg_destruction_game$coefficients[3],
               diff_african_placebo_estimate, diff_western_placebo_estimate,
               african_countries$coefficients[2], african_countries$coefficients[3],
               western_countries$coefficients[2], western_countries$coefficients[3]),
  lower = c(reg_sd$conf.low[2], reg_sd$conf.low[3],
            reg_ss_rights_index$conf.low[2], reg_ss_rights_index$conf.low[3],
            reg_ss_morality_index$conf.low[2], reg_ss_morality_index$conf.low[3],
            reg_message$conf.low[2], reg_message$conf.low[3],
            reg_destruction_game$conf.low[2], reg_destruction_game$conf.low[3],
            diff_african_placebo_estimate - qnorm(1-.05/2)*diff_african_placebo_se, 
            diff_western_placebo_estimate - qnorm(1-.05/2)*diff_western_placebo_se,
            african_countries$conf.low[2], african_countries$conf.low[3], 
            western_countries$conf.low[2], western_countries$conf.low[3]),
  upper = c(reg_sd$conf.high[2], reg_sd$conf.high[3],
            reg_ss_rights_index$conf.high[2], reg_ss_rights_index$conf.high[3],
            reg_ss_morality_index$conf.high[2], reg_ss_morality_index$conf.high[3],
            reg_message$conf.high[2], reg_message$conf.high[3],
            reg_destruction_game$conf.high[2], reg_destruction_game$conf.high[3],
            diff_african_placebo_estimate + qnorm(1-.05/2)*diff_african_placebo_se, 
            diff_western_placebo_estimate + qnorm(1-.05/2)*diff_western_placebo_se,
            african_countries$conf.high[2], african_countries$conf.high[3], 
            western_countries$conf.high[2], western_countries$conf.high[3])
)

# reorder outcomes for the visualization
main_viz$outcome <- factor(main_viz$outcome, 
                           levels = c("Western Country Perceptions",
                                      "African Country Perceptions",
                                      "Double List",
                                      "Contribution Game",
                                      "MP Message",
                                      "Morality Index",
                                      "Rights Index",
                                      "Social Distance"))

# make visualization
ggplot(main_viz, aes(color = type)) + 
  geom_vline(aes(xintercept = 0), lty = "dashed") +
  geom_errorbarh(aes(xmin = lower, xmax = upper, y = outcome),
                 position = position_dodge(width = -0.5),
                 height = 0, size = 1) +
  geom_point(aes(x = estimate, y = outcome, shape = type),
             position = position_dodge(width = -0.5), size = 2.5) + 
  scale_color_manual(values = c("black", "darkgray")) +
  labs(x = "\n Estimate",
       y = "Outcome \n") +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1),
        axis.text=element_text(size=12),
        axis.title=element_text(size=14, face="bold"),
        legend.title = element_blank()) +
  ggsave(filename = "effects_plot.pdf", 
         device = "pdf", width = 8, height = 6, units = "in")

